%% This m-file is to display the reflectivity data
close all
clear all
clc
RefTxtFile={
   'sph01_112_ref.txt';   %   #1 sph01.dat SPM+ 100 mM KCl; pH 3;
    'sph01_141_ref.txt';   %  #2 sph01.dat SPM+ 100 mM KCl; pH 4;
    'sph01_177_ref.txt';   %  #3 sph01.dat SPM/Fe+ 100 mM KCl; pH 3;
    'sph01_251_ref.txt';   %  #4 sph01.dat SPM/Fe+ 100 mM KCl; pH 7.08;
    'sph01_321_ref.txt';   %  #5 sph01.dat SPM/Fe+ 100 mM KCl; pH 4;
    'sph02_27_ref.txt';   %   #6 sph02.dat SPM/Fe+ 100 mM KCl; pH 2;
    'sph02_92_ref.txt';   %   #7 sph02.dat SPM/Fe+ 100 mM KCl; pH 5;
    'sph02_152_ref.txt';   %   #8 sph02.dat SPM/Fe+ 100 mM KCl; pH 5;
    'sph02_215_ref.txt';   %   #9 sph02.dat SPM/Fe+ 100 mM KCl; pH 3;
    'sph02_274_ref.txt';   %   #10 sph02.dat SPM/Fe+ 100 mM KCl; pH neutral 
    'sph02_333_ref.txt';   %   #11 sph02.dat SPM/Fe+ 100 mM KCl; pH neutral 
    'sph02_341_ref.txt';   %   #12 sph02.dat SPM/Fe+ 100 mM KCl; pH neutral 
    }

OptimParamSave={
    'sph01_112_ref_optim.txt';   %  #1  sph01.dat SPM    + 100 mM KCl; pH 3
    'sph01_141_ref_optim.txt';   %  #2  sph01.dat SPM    + 100 mM KCl; pH 4
    'sph01_177_ref_optim.txt';   %  #3  sph01.dat SPM/Fe + 100 mM KCl; pH 3
    'sph01_251_ref_optim.txt';   %  #4  sph01.dat SPM/Fe + 100 mM KCl; pH 7.08
    'sph01_321_ref_optim.txt';   %  #5  sph01.dat SPM/Fe + 100 mM KCl; pH 4
    'sph02_27_ref_optim.txt';    %  #6  sph02.dat SPM/Fe + 100 mM KCl; pH 2
    'sph02_92_ref_optim.txt';    %  #7  sph02.dat SPM/Fe + 100 mM KCl; pH 5
    'sph02_152_ref_optim.txt';   %   #8  sph02.dat SPM/Fe + 100 mM KCl; pH 5
    'sph02_215_ref_optim.txt';   %   #9  sph02.dat SPM/Fe + 100 mM KCl; pH 3
    'sph02_274_ref_optim.txt';   %   #10 *  sph02.dat SPM/Fe + 100 mM KCl; pH neutral
    'sph02_333_ref_optim.txt';   %   #11 * sph02.dat SPM/Fe + 100 mM KCl; pH neutral
    'sph02_341_ref_optim.txt';   %   #12  sph02.dat SPM/Fe + 100 mM KCl; pH neutra
    }

COLOR={'k','b','r',[0 0.4 0],'m',[0.4 0.4 0],[0 0.4 0.6],'b','r'}
Marker={'o','s','v','^','<','>','d'};

C_length=10;
C=zeros(1,C_length);
C(1)=150; % 100 slices
C(6)=0;
C(5)=6.8288; % 2*pi/lambda
%%

PlotRange=[1 6 3 5  ] 
% PlotRange=[1 7 4]
% 1-> pH 3 no Fe
% 6-> pH 2 
% 3-> pH 3
% 5-> pH 4
% 7-> pH 5
% 4-> pH 7

close all
figure
subplot(1,2,1)
hold on
counter=1;
C(2)=0;
for k=PlotRange
    qfit=(0:0.0001:0.55)';
    P=dlmread(OptimParamSave{k});
    dq=P(1);
    I0=P(3);
  
    rfresnel_fit=FresnelReflection(qfit-dq,0.0217,0,1.54);
    rof_fit=fitref_ED_3box(P,C,qfit)./rfresnel_fit/I0;

    data=dlmread(RefTxtFile{k});
    q=data(:,1);
    ref=data(:,2);
    sref=data(:,3);
   
    rfresnel=FresnelReflection(q-dq,0.0217,0,1.54);
    rof=ref./rfresnel/I0;

    plot(qfit,rof_fit,'-','LineWidth',3,'Color',COLOR{counter});
    plot(q,rof,'Marker',Marker{counter},'MarkerFaceColor','w',...
        'MarkerSize',12,'Color',COLOR{counter},...
        'LineWidth',3,'linestyle','none');
    counter=counter+1;
end
hold off
set(gca,'yscale','log');
set(gca,'LineWidth',3);
set(gca,'xlim',[0 0.56]);
%set(gca,'xtick',[0:0.05:0.25]);
set(gca,'xtick',[0:0.1:0.5]);
set(gca,'ylim',[0.5e-2,1e1]);
% set(gca,'xtick',[0:0.1:0.5]);
% set(gca,'ylim',[0.5e-1,1e1]);
set(gca,'box','on');
set(gca,'FontSize',30);


subplot(1,2,2);
hold on
counter = 1;
for k = PlotRange
    z = (-200:1:50)';  % z range
    C(2) = 1;
    P = dlmread(OptimParamSave{k});
    Pnew = P;
    Pnew(7) = P(4);
    Pnew(10) = P(4);
    Pnew(13) = P(4);
    rho = fitref_ED_3box(P, C, z);  % Calculate ED profile
    rhonew = fitref_ED_3box(Pnew, C, z);  % Calculate ED profile
    A = sum(rho-rhonew)

    % Find the point where the ED profile first exceeds 0.1
    I = length(rho);
    for kk = length(rho):-1:1
        if rho(kk) > 0.1
            I = kk; break;
        end
    end
    
    znew = z - z(I);  % Adjust z to match this shift
    plot(znew, rho, '-', 'LineWidth', 3, 'Color', COLOR{counter});
    plot(znew, rhonew, '-.', 'LineWidth', 3, 'Color', COLOR{counter});
    
    % If the first dataset, also calculate bulk and surface rho
    if k == 1
        P_bulk = P;
        P_bulk(7:3:13) = 0;
        rho_bulk = fitref_ED_3box(P_bulk, C, z);
        plot(znew, rho_bulk, '--', 'linewidth', 3, 'Color', COLOR{counter});
       
        rho_surf = rho - rho_bulk;
        plot(znew, rho_surf, '--', 'linewidth', 3, 'Color', COLOR{counter});
    end
    
    % Subtract the constant value 0.334 from rho before integration
    rho_adjusted = rho - 0.334;
    
    % Limit the range for integration from -80 to 0
    valid_range = (znew >= -80) & (znew <= 0);  % Only consider z values from -80 to 0
    z_valid = znew(valid_range);  % Corresponding z values
    rho_valid = rho_adjusted(valid_range);  % Corresponding adjusted rho values
    
    % Check if the range is valid before integrating
    if ~isempty(z_valid) && ~isempty(rho_valid)
        % Interpolate the data for better integration with integral
        rho_interp = @(zq) interp1(z_valid, rho_valid, zq, 'linear', 'extrap');
        
        % Integrate using integral function over the range -80 to 0
        area_under_curve = integral(rho_interp, -80, -10);
        
        % Display the adjusted area in the command window
        fprintf('Adjusted area for dataset #%d: %.4f\n', k, area_under_curve);
    else
        fprintf('No valid data for integration in dataset #%d\n', k);
    end
    % % Save rho and znew to a text file
    % filename = sprintf('rho_z_dataset_%d.txt', k);  % Generate a filename for each dataset
    % fileID = fopen(filename, 'w');  % Open the file for writing
    % 
    % fprintf(fileID, 'znew\t\trho\n');  % Write header
    % 
    % for i = 1:length(znew)
    %     fprintf(fileID, '%f\t%f\n', znew(i), rho(i));  % Write znew and rho values
    % end
    % 
    % fclose(fileID);  % Close the file
    
    counter = counter + 1;
end
% Add vertical line at x = -10
xline(-10, '--r', 'LineWidth', 2);  % '--r' adds a red dashed line
hold off

set(gca, 'LineWidth', 3);
set(gca, 'xlim', [-85, 25]);
set(gca, 'xtick', [-80:20:20]);
set(gca, 'ylim', [-0.1, 0.71]);
set(gca, 'ytick', [0:0.1:0.7]);
set(gca, 'box', 'on');
set(gca, 'FontSize', 30);


